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The observation of massive black hole binaries (MBHBs) with Pulsar Timing Arrays (PTAs) is 
one of the goals of gravitational wave astronomy in the coming years. Massive ( ^ 10** Mq) and 
low-redshift ( ^ 1.5) sources are expected to be individually resolved by up-coming PTAs, and our 
ability to use them as astrophysical probes will depend on the accuracy with which their parameters 
can be measured. In this paper we estimate the precision of such measurements using the Fisher- 
information- matrix formalism. For this initial study we restrict to "monochromatic" sources, i.e. 
binaries whose frequency evolution is negligible during the expected « 10 yr observation time, 
which represent the bulk of the observable population based on current astrophysical predictions. 
In this approximation, the system is described by seven parameters and we determine their expected 
statistical errors as a function of the number of pulsars in the array, the array sky coverage, and the 
signal-to-noise ratio (SNR) of the signal. At fixed SNR (regardless of the number of pulsars in the 
PTA), the gravitational wave astronomy capability of a PTA is achieved with ~ 20 pulsars; adding 
more pulsars (up to 1000) to the array reduces the source error-box in the sky AQ by a factor 
~ 5 and has negligible consequences on the statistical errors on the other parameters, because the 
correlations amongst parameters are already removed to a large extend. If one folds in the increase 
of coherent SNR proportional to the square root of the number of pulsars. Ail improves as 1/SNR^ 
and the other parameters as 1/SNR. For a fiducial PTA of 100 pulsars uniformly distributed in the 
sky and a coherent SNR — 10, we find Af2 ~ 40deg^, a fractional error on the signal amplitude 
of ~ 30% (which constraints only very poorly the chirp mass - luminosity distance combination 
M^^'-^/Dl), and the source inclination and polarization angles are recovered at the ~ 0.3 rad level. 
The ongoing Parkes PTA is particularly sensitive to systems located in the southern hemisphere, 
where at SNR = 10 the source position can be determined with AQ, ~ lOdeg^, but has poorer (by 
an order or magnitude) performance for sources in the northern hemisphere. 

PACS numbers; 



I. INTRODUCTION 



Pulsar Timing Arrays (PTAs), such as the Parkes 
PTA (PPTA) [H, the European PTA (EPTA) 0, 
Nanograv the International Pulsar Timing Array 
(IPTA) project [4L, and in the future the Square Kilome- 
tre Array (SKA) ^ provide a unique means to study the 
population of massive black hole (MBH) binary systems 
with masses above ~ 10'' Mq by monitoring stable radio 
pulsars: in fact, gravitational waves (GWs) generated by 
MBH binaries (MBHBs) affect the propagation of elec- 
tromagnetic signals and leave a distinct signature on the 
time of arrival of the radio pulses . MBH formation 
and evolution scenarios |10l-[l3l| predict the existence of 
a large number of MBHBs. Whereas the high redshift, 
low(er) mass systems will be targeted by the plan ned 
Laser Interferometer Space Antenna {LISA [31) |l5l - [l9j . 
massive and lower redshift {z ^ 2) binaries radiating in 
the (gravitational) frequency range ~ 10~^ Hz — 10~^ Hz 
will be directly accessible to PTAs. These systems im- 
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print a typical signature on the time-of-arrival of radio- 
pulses at a level of w 1 — 100 ns [2^, which is comparable 
with the timing stability of several pulsars [2l| , with more 
expected to be discovered and monitored in the future. 
PTAs therefore provide a direct observational window 
onto the MBH binary population, and can contribute to 
address a number of astrophysical open issues, such as 
the shape of the bright end of the MBH mass function, 
the nature of the MBH-bulge relation at high masses, and 
the dynamical evolution at sub-parsec scales of the most 
massive binaries in the Universe (particularly relevant to 
the so-called "final parsec problem" (22j). 

Gravitational radiation from the cosmic population of 
MBHBs produces two classes of signals in PTA data: (i) 
a stochastic GW background generated by the incoherent 
superposition of radiation from the whole MBHB popu- 
lation [23l - [28| and (ii) individually resolvable, determin- 
istic signals produced by single sources that are suffi- 
ciently massive and/or close so that the gravitational sig- 
nal stands above the root-mean-square (rms) level of the 
background In (SVV, hereafter) we explored 

a comprehensive range of MBH population models and 
found that, assuming a simple order-of-magnitude crite- 
rion to estimate whether sources are resolvable above the 
background level « 1-to-lO individual MBHBs could be 
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observed by future PTAs surveys. The observation of 
GWs from individual systems would open a new avenue 
for a direct census of the properties of MBHBs, offer- 
ing invaluable new information about galaxy formation 
scenarios. The observation of systems at this stage along 
their merger path would also provide key insights into the 
understanding of the interaction between MBHBs and 
the stellar/gaseous environment and how these in- 
teractions affect the black hole-bulge correlations during 
the merger process. If an electro-magnetic counterpart 
of a MBHB identified with PTAs was to be found, such a 
system could offer a unique laboratory for both accretion 
physics (on small scales) and the interplay between black 
holes and their host galaxies (on large scales). 

The prospects of achieving these scientific goals raise 
the question of what astrophysical information could be 
extracted from PTA data and the need to quantify the 
typical statistical errors that will affect the measure- 
ments, their dependence on the total number and spa- 
tial distribution of pulsars in the array (which affects the 
surveys observational strategies), and the consequences 
for multi-band observations. In this paper we estimate 
the statistical errors that affect the measurements of the 
source parameters focusing on MBHBs with no spins, 
in circular orbits, that are sufRciently far from coales- 
cence so that gravitational radiation can be approxi- 
mated as producing a signal with negligible frequency 
drift during the course of the observation time, T w 10 
yr ("monochromatic" signal). This is the class of sig- 
nals that in SVV we estimated to produce the bulk of 
the observational sample. The extension to eccentric bi- 
naries and systems with observable frequency derivative 
is deferred to a future work. GWs from monochromatic 
circular binaries constituted by non-spinning MBHs are 
described by seven independent parameters. We compute 
the expected statistical errors on the source parameters 
by evaluating the variance-covariance matrix - the in- 
verse of the Fisher information matrix - of the observable 
parameters. The diagonal elements of such a matrix pro- 
vide a robust lower limit to the statistical uncertainties 
(the so-called Cramer- Rao bound [10, [3l| ) , which in the 
limit of high signal-to-noise ratio (SNR) tend to the ac- 
tual statistical errors. Depending on the actual structure 
of the signal likelihood function and the SNR this could 
underestimate the actual errors, see e.g. [32| - [3^ for a dis- 
cussion in the context of GW observations. Nonetheless, 
this analysis serves as an important benchmark and can 
then be refined by carrying out actual analyses on mock 
data sets and by estimating the full set of (marginalised) 
posterior density functions of the parameters. The main 
results of the paper can be summarised as follows: 

• At least three (not co-aligned) pulsars in the PTA 
are necessary to fully resolve the source parameters; 

• The statistical errors on the source parameters, at 
fixed SNR, decrease as the number of pulsars in the 
array increases. The typical accuracy greatly im- 
proves by adding pulsars up to ~ 20; for larger ar- 



rays, the actual gain become progressively smaller 
because the pulsars "fill the sky" and the effective- 
ness of further triangulation saturates. In partic- 
ular, for a fiducial case of an array of 100 pulsars 
randomly and uniformly distributed in the sky with 
optimal coherent SNR = 10 ~ which may be appro- 
priate for the SKA - we find a typical GW source 
error box in the sky w 40 dcg^ and a fractional 
amplitude error of w 30%. The inclination and po- 
larization angles can be determined within an er- 
ror of ~ 0.3 rad, and the (constant) frequency is 
determined to sub-frequency resolution bin accu- 
racy. These results are independent on the source 
gravitational- wave frequency. 

• When an anisotropic distribution of pulsars is con- 
sidered, the typical source sky location accuracy 
improves linearly with the array sky coverage. The 
statistical errors on all the other parameters are 
essentially insensitive to the PTA sky coverage, as 
long as it covers more than 1 srad. 

• The ongoing Parkes PTA aims at monitoring 20 
pulsars with a 100 ns timing noise; the targeted 
pulsars are mainly located in the southern sky. A 
GW source in that part of the sky could be local- 
ized down to a precision of < 10 dcg^ at SNR= 10, 
whereas in the northern hemisphere, the lack of 
monitored pulsars limits the error box to }t 200 
deg^. The median of the Parkes PTA angular res- 
olution is « 130 (SNR/10)-2 deg2. 

The paper is organised as follows. In Section II we 
describe the GW signal relevant to PTA and we intro- 
duce the quantities that come into play in the parameter 
estimation problem. A review of the Fisher information 
matrix technique and its application to the PTA case are 
provided in Section HI. Section IV is devoted to the de- 
tailed presentation of the results, and in Section V we 
summarize the main findings of this study and point to 
future work. Unless otherwise specified, throughout the 
paper we use geometric units G = c = I. 

II. THE SIGNAL 

Observations of GWs using PTAs exploit the regular- 
ity of the time of arrival of radio pulses from pulsars. 
Gravitational radiation affects the arrival time of the 
electromagnetic signal by perturbing the null geodesies 
of photons traveling from a pulsar to the Earth. This 
was realized over thirty years ago and the num- 

ber and timing stability of radio pulsars known today 
and expected to be monitored with future surveys [l|- 
H, Q make ensembles of pulsars - PTAs - "cosmic de- 
tectors" of gravitational radiation in the frequency range 
~ 10~® Hz — 10~^ Hz. Here we review the signal produced 
by a GW source in PTA observations. Let us consider 
a GW metric perturbation hab{t) in the transverse and 
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traceless gauge (TT) described by the two independent 
(and time-dependent) polarisation amplitudes h^[t) and 
hy,{t) that carry the information about the GW source. 
Let us also indicate with Cl the unit vector that identifies 
the direction of GW propagation (conversely, the direc- 
tion to the GW source position in the sky is —Cl). The 
metric perturbation can therefore be written as: 

hab{t,(l) ^ e+^^{9,)h+{t,(l) + e'^^{(l)hy,{t,(l), (1) 

where e^{,(f2) {A — + , x) are the polarisation tensors, 
that are uniquely defined once one specifies the wave 
principal axes described by the unit vectors m and h 
as. 



^tbi^) = "^alTib - flaflb , 
'^abi^) = IT^'aflb + flarhb . 



(2a) 
(2b) 



Let us now consider a pulsar emitting radio pulses with a 
frequency i^q. Radio waves propagate along the direction 
described by the unit vector p, and in the background hab 
the frequency of the pulse is affected. For an observer at 
Earth (or at the Solar System Barycentre), the frequency 
is shifted according to the characteristic two-pulse func- 
tion d 



Ahabit;Cl). 



(3) 



Here i/(t) is the received frequency (say, at the Solar Sys- 
tem Barycentre), and 



Ahab{t) = hab{tp,n) ^ hab{t,Cl) 



(4) 



is the difference between the metric perturbation at the 
pulsar - with spacetime coordinates {tp, Xp) - and at the 
receiver - with spacetime coordinates (t,x). The quan- 
tity that is actually observed is the time-residual r{t), 
which is simply the time integral of Eq. ([3]) , 



r{t) = / dt'z{t\n) . (5) 

"'0 



We can re- write Eq. ([3]) in the form 

z{t, ^) = Y^ F'^{n)AhAit; Cl) , 



where 



(6) 



(7) 



^ 1 +p''ila 



is the "antenna beam pattern", see Eqs ([T|), ((2a|) 
and (j2b[) ): here we use the Einstein summation con- 
vention for repeated indeces. Using the definitions (|2ap 
and (j2bp for the wave polarisation tensors, it is simple 



to show that the antenna beam patterns depend on the 
three direction cosines rh ■ p, n ■ p and Cl ■ p: 



F^(J7) 



1 {m ■ p)^ — [h ■ pY 



2 l + n-p 
(m ■ p) [h ■ p) 



l+Vl-p 



(8b) 



Let us now consider a reference frame {x,y, z) fixed to 
the Solar System Barycentre. The source location in the 
sky is defined by the usual polar angles {9, (j)). The unit 
vectors that define the wave principal axes are given by 
(cf. Eqs. (B4) and (B5) in Appendix B of [111; here we 
adopt the same convention used in high-frequency laser 
interferometric observations) 

TO = (sin (j) cos "0 — sin tJj cos cos 9)x 
— (cos 4> COS tJj + sin tJj sin (f> cos 9)y 
+ {sin tp sin 9) z , (9a) 
n ~ (— sin (j) sin tjj — cos if) cos (j) cos 9)x 
+ (cos (j) sin ip — cos iji sin (j) cos 9)y 
+ {costpsin9)z , (9b) 

where i, y and z are the unit vectors along the axis of the 
reference frame, x, y, and z, respectively. The angle tp is 
the wave polarisation angle, defined as the angle counter- 
clockwise about the direction of propagation from the line 
of nodes to the axis described by m. The wave propagates 
in the direction = to x n, which is explicitly given by: 

Vl = — (sin 6' cos 0) x — (sin 6* sin 0) y — cos 6*2 . (10) 

Analogously, the unit vector 

Pa ~ (sin 9a cos (f>a) x + {sin 9a sin (f>a) y + cos 9aZ (11) 

identifies the position in the sky of the a-th pulsar using 
the polar angles {9a, (pa)- 

We will now derive the expression of the PTA signal, 
Eq. [21 produced by a circular, non-precessing binary sys- 
tem of MBHs emitting almost monochromatic radiation, 
i.e with negligible frequency drift during the observa- 
tion time, r « 10 yr. The results are presented in Sec- 
tion IIIBI In the next sub-section, we firstly justify the 
astrophysical assumptions. 



A. Astrophysical assumptions 

Let us justify (and discuss the limitations of) the 
assumptions that we have made on the nature of the 
sources, that lead us to consider circular, non-precessing 
binary systems generating quasi-monochromatic radia- 
tion, before providing the result in Eq. ([31]) . We derive 
general expressions for the phase evolution displacement 
introduced by the frequency drift and by the eccentric- 
ity induced periastron precession, and the change in the 
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orbital angular momentum direction caused by the spin- 
orbit coupling induced precession. The size of each of 
these effects is then evaluated by considering a realis- 
tic (within our current astrophysical understanding) se- 
lected population of resolvable MBHBs taken from SVV. 
Throughout the paper we will consider binary systems 
with masses toi and TO2 ("t-2 "tii), and chirp mass 
Ai = rr?-^^ I {mx -^-m^)^!'^ , emitting a GW frequency 
/. We also define M — m\ -I- 7712, /i = m\milm and 
q = m2/mi, the total mass, the reduced mass and the 
mass ratio, respectively. Our notation is such that all the 
quantities are the observed (redshifted) ones, such that 
e.g. the intrinsic (rest-frame) mass of the primary MBH 
is mi^r = + z) and the rest frame GW frequency 

is fr = /(I + z). We normalize all the results to 



Ms.5 

/50 = 



M 



109 Mq 
M 



108-5 Mq ' 
/ 



50 nHz ' 
T 



10 yr 



which are the typical values for individually resolvable 
sources found in SVV, and the typical observation times- 
pan. 



velocity of the binary is: 

V ^ (7r/il/)2/3 



1.73 X 10-2Af2/3/2/3 



'50 



(14) 



Stated in different terms, the systems will be far from 
plunge, as the time to coalescence for a binary radiating 
at frequency / is (at the leading Newtonian quadrupole 
order and for a circular orbit system) 



ieoal=^4x 10^7W8.^'^/. 



5/3 



yr- 



(15) 



As a consequence the frequency evolution during the ob- 
servation time is going to be small, and can be neglected. 
In fact, it is simple to estimate the total frequency shift 
of radiation over the observation period 



A/ « /T « 0.05 Ml'l fll'^ Tio nHz 



(16) 



which is negligible with respect to the frequency reso- 
lution bin « 3r^Q^ nHz; correspondingly, the additional 
phase contribution 



A<I> « tt/T^ « 0.04A^^/5' fll'^ rad 



(17) 



is much smaller than 1 rad. Eqs. and p7)) clearly 
show that it is more than legitimate in this initial study to 
ignore any frequency derivative, and treat gravitational 
radiation as monochromatic over the observational pe- 
riod. 



1. Gravitational wave frequency evolution 

A binary with the properties defined above, evolves 
due to radiation reaction through an adiabatic in-spiral 
phase, with GW frequency f(t) changing at a rate (at 
the leading Newtonian order) 

^^^^8/3^5/3^11/3. (12) 

The in-spiral phase terminates at the last stable orbit 
(LSO), that for a Schwarzschild black hole in circular 
orbit corresponds to the frequency 

/lso =4.4x lO-^Mg-i Hz. (13) 

The observational window of PTAs is set at low fre- 
quency by the overall duration of the monitoring of pul- 
sars T « 10 yr, and at high frequency by the cadence of 
the observation, « 1 week: the PTA observational win- 
dow is therefore in the range ^ 10^^ — 10^^ Hz. In SVV 
we explored the physical properties of MBHBs that arc 
likely to be observed in this frequency range: PTAs will 
resolve binaries with mi, 2 ~ lO^M© and in the frequency 
range « 10"* — 10~^ Hz. In this mass-frequency region, 
PTAs will observe the in-spiral portion of the coalescence 
of a binary system and one can ignore post-Newtonian 
corrections to the amplitude and phase evolution, as the 



2. Spin effects 

Wc now justify our assumption of neglecting the spins 
in the modelling of the waveform. From an astrophysi- 
cal point of view, very little precise information about 
the spin of MBHs can be extracted directly from ob- 
servations. However, several theoretical arguments sup- 
port the existence of a population of rapidly spinning 
MBHs. If coherent accretion from a thin disk [3g| is the 
dominant growth mechanism, then MBH spin-up is in- 
evitable [37|; jet production in active galactic nuclei is 
best explained by the presence of rapidly spinning MBHs 
[38| ; in the hierarchical formation context, though MBHB 
mergers tend to spin down the remnant [39| . detailed 
growth models that take into account both mergers and 
accretion lead to populations of rapidly spinning MBHs 
[40I l4l| . Spins have two main effects on the gravitational 
waveforms emitted during the in-spiral: (i) they affect 
the phase evolution [S^l, and (ii) they cause the orbital 
plane to precess through spin-orbit and spin-spin cou- 
pling [431 . Il^ l . The effect of the spins on the phase evolu- 
tion is completely negligible for the astrophysical systems 
observable by PTAs: the additional phase contribution 
enters at the lowest order at the post^ 5_]vj^e.^^,^Qnia,n order, 
that is proportional to v'^ , and we have already shown 
that ^ 1, see Eq. ((Ti)) . Precession would provide a 
characteristic imprint on the signal through amplitude 
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and phase modulations produced by the orbital plane 
precession, and as a consequence the time-dependent po- 
larisation of the waves as observed by a PTA. It is fairly 
simple to quantify the change of the orientation of the or- 
bital angular momentum unit vector L during a typical 
observation. The rate of change of the precession angle 
is at the leading order: 



~dt 



= 2 



37Tl2 \ L + S 

2mi J a-^ 



(18) 



where L = ^ a^^M is the magnitude of the orbital an- 
gular momentum and 5* is the total intrinsic spin of the 
black holes. As long as fJ,/M ^ w/c, we have that L ^ S. 
This is always the case for resolvable MBHBs; we find in- 
deed that these systems are in general characterised by 
q > 0.1 (therefore fi/M > 0.1), while from Eq. ([M]) we 
know that in general v/c ^ 0.01 . In this case, from 
Eq. ^TE\\ one obtains 



2^5/3 1 



0.8 1 + 



4m 1 
3m2 
4toi 



which is independent of 5*. The effect is maximum for 
equal mass binaries, mi = m2, [ijM = 0.25; in this case 



Ac 



0.3 rad. It is therefore clear that in general spins 



will not play an important role, and wc will neglect their 
effect in the modeling of signals at the PTA output. It 
is however interesting to notice that for a IQ^Mq binary 
system observed for 10 years at « 10"'' Hz, which is con- 
sistent with astrophysical expectations (see SVV) the ori- 
entation of the orbital angular momentum would change 
by Aofp « 1 rad. The Square-Kilometre- Array has there- 
fore a concrete chance of detecting this signature, and to 
provide direct insights onto MBH spins. 



3. Eccentricity of the binary 

Let us finally consider the assumption of circular or- 
bits, and the possible effects of neglecting eccentricity in 
the analysis. The presence of a residual eccentricity at or- 
bital separations corresponding to the PTA observational 
window has two consequences on the observed signal: (i) 
the power of radiation is not confined to the harmonic 
at twice the orbital frequency but is spread on the (in 
principle infinite) set of harmonics at integer multiples 
of the inverse of the orbital period, and (ii) the source 
periapsc precesses in the plane of the orbit at a rate 

~ 3.9 X 10-3 (1 - e^) Mg /g rad s'^ (20) 

which introduces additional modulations in phase and (as 
a consequence) amplitude in the signal recorded at the 



Earth. In Eq. (|20)) 7(t) is the angle of the periapse mea- 
sured with respect to a fixed frame attached to the source. 
We now briefiy consider the two effects in turn. The pres- 
ence of eccentricity "splits" each polarisation amplitude 
h^[t) and h^it) into harmonics according to (see e.g. 
Eqs. (5-6) in Ref. [i^ and references therein): 

h+{t) ^ yl|-(l + cos2i)u„(e)cos ^ $(t) -I- 27(t) 



-(1 + cos^ t)'^n(e) cos -<I>(i)-27(t) 



where 



-I- sin^ iWn{e) cos 
2 A cos t|w„(e) su 

{ft 1 "1 

+«„(e)sin([-$(t)-27(i)J)| 



-$(0 + 27(0 



(21) 



(22) 



m=2^ f fit')dt' , (23) 

is the GW phase and f{t) the instantaneous GW fre- 
quency corresponding to twice the inverse of the or- 
bital period. The source inclination angle l is defined 
as cosi = —fl^'La, where L is the unit vector that de- 
scribes the orientation of the source orbital plane, and 
the amplitude coefficients Mn(e), Wn(e), and w„(e) are 
linear combinations of the Bcssel functions of the first 
kind J„(ne), Jn±i{ne) and Jn±2{'ne). For an astrophysi- 
cally plausible range of eccentricities e 0.3 - see Fig. [1] 
and the discussion below - |Mn(e)| ^ kn(e)| , |w„(e)| and 
most of the power will still be confined into the n = 2 
harmonic at twice the orbital frequency, see e.g. Fig. 3 
of Ref. [IHl- On the other hand, the change of the peri- 
apse position even for low eccentricity values may intro- 
duce significant phase shifts over coherent observations 
lasting several years. In fact the phase of the recorded 
signal is shifted by an additional contribution 27(0. This 
means that the actual frequency of the observed signal 
recorded at the instrument corresponds to f{t) + 'j/TT and 
differs by a measurable amount from f{t). Nonetheless, 
one can still model the radiation observed at the PTA 
output as monochromatic, as long as th periapse preces- 
sion term 7/71 introduces a phase shift quadratic in 
time that is ^ 1 rad, which is equivalent to the condi- 
tion that we have imposed on the change of the phase 
produced by the frequency shift induced by radiation re- 
action, see Eqs. ^ and From Eq. ^ and ([12]), 
this condition yields: 



A$^ 



dt^ (1 — e^) 



d'1^2 



2 X 10^^(1 -e^) 



(24) 



We therefore see that the effect of the eccentricity will be 
in general negligible. 
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A4> [rad] Aa [rad] 

lO-nO-nO-^O.l 1 10-3 iQ-z 0.1 1 




10-3 10-2 0.1 
eccentricity 



1 10-s 10-" 10-3 10-2 0.1 
A<J>^ [rad] 



FIG. 1: Testing the circular monochromatic non-spinning bi- 
nary approximation. Upper left panel: distribution of the 
phase displacement A$ introduced by the frequency drift of 
the binaries. Upper right panel: change in the orbital an- 
gular momentum direction Aop introduced by the spin-orbit 
coupling. Lower left panel: eccentricity distribution of the 
systems. Lower right panel: distribution of phase displace- 
ment A"I>^ induced by relativistic periastron precession due 
to non-zero eccentricity of the binaries. The distributions are 
constructed considering all the resolvable MBHBs with resid- 
uals > Ins (solid lines), 10ns (long-dashed lines) and 100ns 
(short-dashed lines), found in 1000 Monte Carlo realizations 
of the Tu-SA models described in SVV, and they are nor- 
malised so that their integrals are unity. 



4- Tests on a massive black hole population 

We can quantify more rigorously whether the assump- 
tion of monochromatic signal at the PTA output is jus- 
tified, by evaluating the distributions of A$, Aap and 
on an astrophysically motivated population of re- 
solvable MBHBs. We consider the Tu-SA MBHB popula- 
tion model discussed in SVV (see Section 2.2 of SVV for 
a detailed description) and we explore the orbital evo- 
lution, including a possible non-zero eccentricity of the 
observable systems. The binaries are assumed to be in 
circular orbit at the moment of pairing and are self con- 
sistently evolved taking into account stellar scattering 
and GW emission 
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We generate 1000 Monte Carlo 
realisations of the entire population of GW signals in 
the PTA band and we collect the individually resolvable 
sources generating coherent timing residuals greater than 
1, 10 and 100 ns, respectively, over 10 years. In Fig. [1] 
we plot the distributions relevant to this analysis. We see 
from the two upper panels, that in general, treating the 
system as "monochromatic" with negligible spin effects 



is a good approximation. If we consider a 1 ns thresh- 
old (solid lines), the phase displacement A$ introduced 
by the frequency drift and the orbital angular momen- 
tum direction change Aa^ due two spin-orbit coupling 
are always < 1 rad, and in ~ 80% of the cases are < 0.1 
rad. The lower left panel of Fig. [1] shows the eccentricity 
distribution of the same sample of individually resolv- 
able sources. Almost all the sources are characterised 
by e ;S 0.1 with a long tail extending down to e ^ lO"'^ 
in the PTA band. The typical periastron precession- 
induced additional phase 2'yT can be larger than 1 rad. 
However, this additional contribution grows linearly with 
time, and, as discussed before, will result in a measured 
frequency which differs from the intrinsic one by a small 
amount j/ir < InHz. The "non-monocromatic" phase 
contribution A$-y that changes quadratically with time 



and is described by Eq. ([24]) is instead plotted in the 
lower right panel of Fig. [TJ Values of A$-y are typically 
of the order lO^'^, completely negligible in the context 
of our analysis. Note that, as a general trend, increasing 
the threshold in the source-induced timing residuals to 10 
and 100 ns, all the effects tend to be suppressed. This is 
because resolvable sources generating larger residuals are 
usually found at lower frequencies, and all the effects have 
a steep dependence on frequency - see Eqs. P?)l . 
and (|24p . This means that none of the effects considered 
above should be an issue for ongoing PTA campaigns, 
which aim to reach a total sensitivity of > 30 ns, but may 
possibly play a role in recovering sources at the level of a 
few ns, which is relevant for the planned SKA. Needless to 
say that a residual eccentricity at the time of pairing may 
result in larger values of e than those shown in Fig. [T] [43], 
causing a significant scatter of the signal power among 
several different harmonics; however, the presence of gas 
may lead to circularization before they reach a frequency 
« 10~^ Hz (see, e.g., HI]). Unfortunately, little is known 
about the eccentricity of subparsec massive binaries, and 
here we tackle the case of circular systems, deferring the 
study of prcccssing eccentric binaries to future work. 



B. Timing residuals 

We have shown that the assumption of circular, 
monochromatic, non-precessing binary is astrophysically 
reasonable, surely for this initial exploratory study. We 
now specify the signal observed at the output, Eq. ([5|), 
in this approximation. The two independent polarisa- 
tion amplitudes generated by a binary system, Eqs. (j2ip 
and (|22|) . can be written as: 



h+{t) = Ag^aii) cos $(t) , 
h^{t) = Agw5(t)sin$(t) , 



where 



^gw(/) 



2^ [^/W] 



(25a) 
(25b) 



(26) 
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Af„ [Hz] 



FIG. 2: Normalized distribution of Afa (see text) for the 
same sample of MBHBs considered in Fig. [TJ assuming ob- 
servations with 100 isotropically distributed pulsars in the 
sky at a distance of 1 kpc. The vertical dotted line marks 
the width of the array's frequency resolution bin Afr = 1/T 



(« 3 X 10"^Hz for T = lOyr) 



is the GW amplitude, D the luminosity distance to the 
GW source, <i>(t) is the GW phase given by Eq. ([23|) . and 
f(t) the instantaneous GW frequency (twice the inverse 
of the orbital period). The two functions 



cos^ L 



a(i) = 1 
6(i) = —2 cos L 



(27a) 
(27b) 



depend on the source inclination angle t, defined in the 
previous Section. 

As described in Section II, Eqs. and (O, the re- 
sponse function of each individual pulsar a consists of 
two terms, namely, the perturbation registered at the 
Earth at the time t of data collection {hab{t,ri)), and 
the perturbation registered at the pulsar at a time t — 
(habit — Ta, f2)), where Tq, is the light-travel-time from the 
pulsar to the Earth given by: 



Ta = La{l + ^- Pa) 
,11 La 



1.1 X 10' 



1 kpc 



(1 + n-pa)s, 



(28) 



where is the distance to the pulsar. We can therefore 
formally write the observed timing residuals, Eq. ([5]) for 
each pulsar a as: 



r„(t) = r(^)(t)-fri^)(t), 



(29) 



where P and E label the "pulsar" and "Earth" contribu- 
tion, respectively. During the time Tq, the frequency of 



the source - although " monochromatic" over the time of 
observation T of several years - changes by 



Afa 



dt" 



dt 



^1 
dt ' 



15 7Wc4V5o^^r„,i nHz, 



(30) 

where r^^i is the pulsar-Earth light-travel-time normal- 
ized to a distance of 1 kpc. The frequency shift Afa 
depends both on the parameters of the source (emission 
frequency and chirp mass) and the properties of the pul- 
sar (distance and sky location with respect to the source). 
We can quantify this effect over an astrophysically plausi- 
ble sample of GW sources by considering the population 
shown in Fig. [T] Let us consider the same set of resolv- 
able sources as above, and assume detection with a PTA 
of 100 pulsars randomly distributed in the sky, but all at 
a distance of 1 kpc. For each source we consider all the 
Afa related to each pulsar and we plot the results in Fig. 
[2l The distribution has a peak around ~ 5 x 10"^ Hz, 
which is ^ 10 times larger than the typical frequency 
resolution bin for an observing time T « 10 yr. This 
means that the signal associated to each pulsar gener- 
ates at the PTA output two monochromatic terms at two 
distinct frequencies. All the "Earth-terms" corresponding 
to each individual pulsar share the same frequency and 
phase. They can therefore be coherently summed among 
the array, building up a distinct monochromatic peak 
which is not going to be affect by the pulsar terms (also 
known as " self- noise" ) which usually happen to be at 
much lower frequencies. The contribution to the Earth 
term from each individual pulsar can be written as 



ri'^Ht) 



R [a F+ (sin <I>(t) - sin $o j 
5F„^(cos$(i) - cos$o) 



with 



(31) 



(32) 



and $(<) given by Eq. (|23p . The Earth timing residuals 
are therefore described by a 7-dimensional vector encod- 
ing all (and only) the parameters of the source: 
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X = {R,9,^,^,i,f, $o}. 



(33) 



Conversely, each individual pulsar term is characterized 
by a different amplitude, frequency and phase, that cru- 
cially depend also on the poorly constrained distance La 
to the pulsar. In order to take advantage of the power 
contained in the pulsar term, one needs to introduce an 
additional parameter for each pulsar in the PTA. As 
a consequence, this turns a 7-parameter reconstruction 
problem into a 7 -|- M parameter problem. More de- 
tails about the PTA response to GWs are given in Ap- 
pendix A. In this paper, we decided to consider simply 
the Earth-term (at the expense of a modest loss in total 
SNR) given by Eq. ([31]) , which is completely specified by 
the 7-parameter vector At present, it is not clear 
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whether it would also be advantageous to include into the 
analysis the pulsar-terms, that require the addition of M 
unknown search parameters. This is an open issue that 
deserves further investigations and will be considered in 
a future paper. 



III. PARAMETER ESTIMATION 

In this section we briefly review the basic theory and 
key equations regarding the estimate of the statistical 
errors that affect the measurements of the source param- 
eters. For a comprehensive discussion of this topic we 
refer the reader to [sot . 

The whole data set collected using a PTA consisting of 
M pulsars can be schematically represented as a vector 



d= {di,d2,...,dM} , 



(34) 



where the data form the monitoring each pulsar (a 
1, . . . , M) are given by 



da{t) = na{t) + ra{t; A) . 



(35) 



In the previous equation ra{t; A), given by Eq. is the 
GW contribution to the timing residuals of the a-th pul- 
sar (the signal) - to simplify notation we have dropped 
(and will do so from now on) the index " E" , but it should 
be understood as we have stressed in the previous sec- 
tion that we will consider only the Earth-term in the 
analysis - and na{t) is the noise that affects the observa- 
tions. For this analysis we make the usual (simplifying) 
assumption that Ha is a zero-mean Gaussian and station- 
ary random process characterised by the one-sided power 
spectral density Sa{f)- 

The inference process in which we are interested in 
this paper is how well one can infer the actual value of 
the unknown parameter vector A, Eq. p3p . based on the 
data d, Eq. p4p . and any prior information on A available 
before the experiment. Within the Bayesian framework, 
see e.g. [4§|, one is therefore interested in deriving the 
posterior probability density function (PDF) p{X\d) of 
the unknown parameter vector given the data set and 
the prior information. Bayes' theorem yields 



piX\d) 



p(A)p(d|A) 
P{d) 



(36) 



where p{d\X) is the likelihood function, p{X) is the prior 
probability density of A, and p{d) is the marginal likeli- 
hood or evidence. In the neighborhood of the maximum- 
likelihood estimate value A, the likelihood function can be 
approximated as a multi-variate Gaussian distribution. 



p{X\d) (X p(A) exp 



iPafcAAaAAb 



(37) 



where AA^ = Xa — Xa and the matrix Tat is the Fisher 
information matrix; here the indexes a, 6 = 1, . . . , 7 la- 
bel the components of A.. Note that we have used Ein- 
stein's summation convention (and we do not distinguish 
between covariant and contravariant indeces). In the 

limit of large SNR, A tends to A, and the inverse of the 
Fisher information matrix provides a lower limit to the 
error covariancc of unbiased estimators of A, the so-called 
Cramer-Rao bound jsij . The variance-covariance matrix 
is simply the inverse of the Fisher information matrix, 
and its elements are 



Cab = 



(38a) 
(38b) 



where — 1 < Cab < +1 (Va, 6) are the correlation coeffi- 
cients. We can therefore interpret cr^ as a way to quan- 
tifying the expected uncertainties on the measurements 
of the source parameters. We refer the reader to [s^ 
and references therein for an in-depth discussion of the 
interpretation of the inverse of the Fisher information 
matrix in the context of assessing the prospect of the es- 
timation of the source parameters for GW observations. 
Here it suffices to point out that MBHBs will likely be 
observed at the detection threshold (see SVV), and the 
results presented in Section ITVl should indeed be regarded 
as lower-limits to the statistical errors that one can ex- 
pect to obtain in real observations, see e.g. [S^-HS- 

One of the parameters that is of particular interest is 
the source sky location, and we will discuss in the next 
Section the ability of PTA to define an error box in the 
sky. Following Ref. [H^l, we define the PTA angular res- 
olution, or source error box as 



An ^ 2TT^{smeAeA(j))^ - (sin6lc^'^)2 ; (39) 

with this definition, the probability for a source to lay 
outside the solid angle AOq is e^-^^o/AO ^ 

We turn now on the actual computation of the Fisher 
information matrix Tab. First of all we note that in ob- 
servations of multiple pulsars in the array, one can safely 
consider the data from different pulsars as independent, 
and the likelihood function of d is therefore 



p{d\x) = n^KiA) 



cx exp 



'-TabAXaAXi, 



(40) 



where the Fisher information matrix that characterises 
the joint observations in the equation above is simply 
given by 



Tab - ^ r^"^^ 



(41) 



r^j'^'' is the Fisher information matrix relevant to the ob- 
servation with the a— th pulsar, and is simply related 
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to the derivatives of the GW signal with respect to the 
unknown parameters integrated over the observation: 



dXa 



drait;X) 



dXh 



(42) 



where the inner product between two functions x(t) and 
y{t) is defined as 



(.|.)^2f ^lMl±im/l./, (43a) 



Snif) 



— / x{t)y{t)dt, 

•JO Jo 



and 



i{f) 



x{t)e 



-2-Kift 



(43b) 



(44) 



is the Fourier Transform of a generic function x{t). The 
second equahty, Eq. (j43b[) is correct only in the case in 
which the noise spectral density is approximately con- 
stant (with value Sq) across the frequency region that 
provides support for the two functions x{f) and y(/). 
Eq. (j43bp is appropriate to compute the scalar product 
for the observation of gravitational radiation from MB- 
HBs whose frequency evolution is negligible during the 
observation time, which is astrophysically justified as we 
have shown in Section ID 

In terms of the inner product (.|.) ~ Eqs. (|43a[) 
and (|43b|) - the optimal SNR at which a signal can be 
observed using a pulsars is 



snr: 



(45) 



and the total coherent SNR produced by timing an array 
of M pulsars is: 



M 



SNR2 ^ SNR^ 



(46) 



IV. RESULTS 

In this section we present and discuss the results of 
our analysis aimed at determining the uncertainties sur- 
rounding the estimates of the GW source parameters. We 
focus in particular on the sky localization of a MBHB, 
which is of particular interest for possible identifica- 
tions of electromagnetic counterparts, including the host 
galaxy and/or galactic nucleus in which the MBHB re- 
sides. For the case of binaries in circular orbit and whose 
gravitational radiation does not produce a measurable 
frequency drift, the mass and distance are degenerate, 
and can not be individually measured: one can only mea- 
sure the combination M.^^^ /Dl. This prevents measure- 
ments of MBHB masses, which would be of great interest. 
On the other hand, the orientation of the orbital angular 



momentum - through measurements of the inclination 
angle l and the polarisation angle ip - can be determined 
(although only with modest accuracy, as we will show be- 
low), which may be useful in determining the geometry 
of the system, if a counterpart is detected. 

The uncertainties on the source parameters depend on 
a number of factors, including the actual MBHB param- 
eters, the SNR, the total number of pulsars and their 
location in the sky with respect to the GW source. It 
is therefore impossible to provide a single figure of merit 
that quantifies of how well PTAs will be able to do GW 
astronomy. One can however derive some general trends 
and scalings, in particular how the results depend on the 
number of pulsars and their distribution in the sky, which 
we call the sky coverage of the array] this is of particu- 
lar importance to design observational campaigns, and 
to explore tradeoffs in the observation strategy. In the 
following subsections, by means of extensive Monte Carlo 
simulations, we study the parameter estimation accuracy 
as a function of the number of pulsars in the array, the 
total SNR of the signal, and on the array sky coverage. 
All our major findings are summarised in Table HI 



A. General behavior 

Before considering the details of the results we dis- 
cuss conceptually the process by which the source pa- 
rameters can be measured. Our discussion is based on 
the assumption that the processing of the data is done 
through a coherent analysis. The frequency of the signal 
is trivially measured, as this is the key parameter that 
needs to be matched in order for a template to remain 
in phase with the signal throughout the observation pe- 
riod. Furthermore, the amplitude of the GW signal de- 
termines the actual SNR, and is measured in a straight- 
forward way. The amplitude R, or equivalently ^gw, see 
Eqs. (pS)) and ([5^ . provides a constraint on the chirp 
mass and distance combination M^^^ /Dl. However, in 
the case of monochromatic signals, these two parame- 
ters of great astrophysical interest can not be measured 
independently. If the frequency derivative / were also 
observable - this case is not considered in this paper, as 
it likely pertains only to a small fraction of detectable 
binaries, see Section [H] and Fig. [5] - then one would be 
able to measure independently both the luminosity dis- 
tance and chirp mass. In fact, from the measurement of 
/ cx 7W5/3jii/3^ that can be evaluated from the phase 
evolution of timing residuals, one can measure the chirp 
mass, which in turn, from the observation of the ampli- 
tude, would yield an estimate of the luminosity distance^. 



^ Wc note that a direct measurement of the chirp mass would be 
possible if one could detect both the Earth- and pulsar-terms, 
provided that the distance to the pulsar was known. In this case 
one has the GW frequency at Earth, the GW frequency at the 
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TABLE I: Typical uncertainties in the measurement of the GW source parameters as a function of the total number of pulsars 
in the array M and their sky coverage AS7pta (the portion of the sky over which the pulsars are uniformly distributed). For 
each PTA configuration we consider 2.5 x 10*-to-1.6 x 10® (depending on the number of pulsars in the array) GW sources with 
random parameters. The GW source location is drawn uniformly in the sky, and the other parameters are drawn uniformly 
over the full range of tp, 4>q and cos t, /o is fixed at 5 x 10~*Hz. In every Monte Carlo realisation, the optimal SNR is equal to 
10. The table reports the median of the statistical errors AA - where A is a generic source parameter - and the 25**^ and 75'i' 
percentile of the distributions obtained from the Monte Carlo samplings. Note that the errors AR/R, Al, Atp, A/, A$o all 
scale as SNR~i, the error ASl scales as SNR^^. 



The remaining parameters, those that determine the ge- 
ometry of the binary ~ the source location in the sky, and 
the orientation of the orbital plane - and the initial phase 
00 can be determined only if the PTA array contains at 
least three (not co-aligned) pulsars. The source location 
in the sky is simply reconstructed through geometrical 
triangulation, because the PTA signal for each pulsar 
encodes the source coordinates in the sky in the relative 
amplitude of the sine and cosine term of the response or, 
equivalently, the overall phase and amplitude of the sinu- 
soidal PTA output signal, see Eqs. ©, dH), dZl) and ([3T1). 
For the reader familiar with GW observations with LISA, 
we highlight a fundamental difference between LISA and 
PTAs in the determination of the source position in the 
sky. With LISA, the error box decreases as the signal fre- 
quency increases (everything else being equal), because 
the source location in the sky is reconstructed (primarily) 
through the location-dependent Doppler effect produced 
by the motion of the instrument during the observation, 
which is proportional to the signal frequency. This is not 
the case for PTAs. where the error-box is independent of 
the GW frequency. It depends however on the number of 
pulsars in the array - as the number of pulsars increases. 



pulsar, and the Earth-pulsar light-travcl-timc, which in turns 
provides a direct measure of /, and as a consequence of the chirp 
mass. 



one has to select with increasingly higher precision the 
actual value of the angular parameters, in order to ensure 
that the same GW signal fits correctly the timing resid- 
uals of all the pulsars - and the location of the pulsars 
in the sky. 

We first consider how the parameter estimation de- 
pends on the total number of pulsars M at fixed SNR. 
We consider a GW source with random parameters and 
we evaluate the inverse of the Fisher information matrix 
as we progressively add pulsars to the array. The pulsars 
are added randomly from a uniform distribution in the 
sky and the noise has the same spectral density for each 
pulsar. We also keep the total coherent SNR fixed, at 
the value SNR = 10. It is clear that in a real observation 
the SNR actually increases approximately as -s/M, and 
therefore depends on the number of pulsars in the array. 
However, by normalising our results to a constant total 
SNR, we are able to disentangle the change in the un- 
certainty on parameter estimation that depends on the 
number of pulsars from the change due simply to the 
SNR. The results are shown in Fig. [3] The main effect 
of adding pulsars in the PTA is to improve the power of 
triangulation and to reduce the correlation between the 
source parameters. At least three pulsars in the array 
are needed to formally resolve all the parameters; how- 
ever, given the strong correlation in particular amongst 
R, L and '0 (which will be discussed later in more de- 
tail) a SNR ~ 100 is needed to locate the source in the 
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FIG. 3: The statistical errors that affect the determination 
of the source location Afi, see Eq. (|39|) (upper panels) and 
the signal amplitude R (lower panels) for four randomly se- 
lected sources (corresponding to the different line styles). We 
increase the number of pulsars in the array fixing a total 
SNR= 10, and we plot the results as a function of the number 
of pulsars M. In the left panels we consider selected edge-on 
(t = it/2) sources, while in the right panel we plot sources 
with intermediate inclination of t = 7r/4. 



sky with an accuracy < 50 deg^ in this case. It is clear 
that the need to maintain phase coherency between the 
timing residuals from several pulsars leads to a steep (by 
orders of magnitude) increase in accuracy from M = 3 
to M w 20 (note that the current Parkcs PTA counts 20 
pulsars). Adding more pulsars to the array reduces the 
uncertainty location region in the sky A51 by a factor of 
« 5 going from 20 to 1000 pulsars, but has almost no im- 
pact on the determination of the other parameters (the 
bottom panels of Fig. [3] show that AR/R is essentially 
constant for M > 20). 

Now that we have explored the effect of the number 
of pulsars alone (at fixed SNR) on the parameter er- 
rors, we can consider the case in which we also let the 
SNR change. We repeat the analysis described above, 
but now the SNR is not kept fixed and we let it vary self- 
consistently as pulsars are added to the array. The results 
plotted as a function of the total coherent SNR are shown 
in Fig. m Once more, we concentrate in particular on 
the measurement of the amplitude R and the error box 
in the sky Afi. For M ^ 1, the error box in the sky and 
the amplitude measurements scale as expected accord- 
ing to An cx SNR-2 and AR/R a SNR"^ (and so do 
all the other parameters not shown here) . However, for 
SNR ^10 the uncertainties departs quite dramatically 
from the scaling above simply due to fact that with only 



FIG. 4: Same as Fig. [3l but here, as we add pulsars to 
the PTA, we consistently take into account the effect on the 
total coherent SNR, and accordingly we plot the results as a 
function of the SNR. In the left panels we plot selected edge- 
on (t = tt/2) sources, while in the right panel we consider 
selected sources with intermediate inclination of i = 7r/4. The 
dotted-dashed thin lines in the upper panels follow the scaling 
AO. oc SNR-^ 



a handful of pulsars in the array the strong correlations 
amongst the parameters degrade the measurements. We 
stress that the results shown here are independent of the 
GW frequency; we directly checked this property by per- 
forming several tests, in which the source's frequency is 
drawn randomly in the range 10^* Hz - 10^^ Hz. 

The source inclination l angle is strongly correlated 
with the signal amplitude R, and the polarisation an- 
gle tp is correlated to both l and <I>o. The results are 
indeed affected by the actual value of the source inclina- 
tion. Left panels in Figs. [3] and |3] refer to four different 
edge-on sources (i.e. l — tt/2 and the radiation is linearly 
polarised). In this case, the parameter have the least cor- 
relation, and AR/R =SNR^^. Right panels in Figs. [3] 
and [3] refer to sources with an "intermediate" inclination 
L = 7r/4; here degeneracies start to play a significant role 
and cause a factor of « 3 degradation in AR/R estima- 
tion (still scaling as SNR~^). Note, however, that the 
sky position accuracy is independent on l (upper panels 
in Figs. [3]and|4|), because the sky coordinates 9 and cj) 
are only weakly correlated to the other source parame- 
ters. We further explore this point by considering the 
behaviour of the correlation coefficients {c^'' and c"^*") 
as a function of l. Fig. \5\ shows the correlation coef- 
ficients and statistical errors in the source's parameters 
for a sample of 1000 individual sources using a PTA with 
AI = 100 and total SNR= 10 as a function of l. For a 
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L [rad] 



FIG. 5: The effect of the source orbital inchnation l on the 
estimate of the signal parameters. Upper panels: The cor- 
relation coefficients c^'' (left) and c'^'^'o (right) as a function 
of L. Middle and bottom panels: the statistical errors in the 
measurement of amplitude R, polarisation angle tp, inclina- 
tion angle and initial phase $o for a fixed PTA coherent SNR 
= 10, making clear the connection between inclination, corre- 
lation (degeneracy) and parameter estimation. Each asterisk 
on the plots is a randomly generated source. 

face-on source (i — 0, tt), both polarizations equally con- 
tribute to the signal, and any polarization angle ip can 
be perfectly 'reproduced' by tuning the source phase $0i 
i.e. the two parameters are completely degenerate and 
cannot be determined. Moving towards edge-on sources, 
progressively change the relative contribution of the two 
polarizations, breaking the degeneracy with the phase. 
Fig. ini shows statistical error distributions for the differ- 
ent parameters over a sample of 25000 sources divided in 
three different t bins. The degradation in the determi- 
nation of R, L and tp moving towards face-on sources is 
clear. Conversely, both 6 and (f> do not have any strongly 
dependent correlation with the other parameters, the es- 
timation of is then independent on the source inclina- 
tion (lower right panel in Fig. IB]). 



B. Isotropic distribution of pulsars 

In this Section we study the parameter estimation for 
a PTA whose pulsars are isotropically distributed in the 
sky, and investigate how the results depend on the num- 
ber M of pulsars in the array and the SNR. Current 
PTAs have pulsars that are far from being isotropically 
located on the celestial sphere - the anisotropic distri- 
bution of pulsars is discussed in the next Section - but 
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FIG. 6: The distributions of the statistical errors of the source 
parameter measurements using a sample of 25000 randomly 
distributed sources (see text for more details), divided in three 
different inclination intervals: t G [0, 7r/6] U [5/6tt, tt] (dotted), 
L G [7r/6,7r/3] U [2/37r, S/Btt] (dashed) and t £ [tt/S, 2/37r] 
(solid). In each panel, the sum of the distribution's integrals 
performed over the three l bins is unity. 



the isotropic case is useful to develop an understanding 
of the key factors that impact on the PTA performances 
for astronomy. It can also be considered representative 
of future PTAs, such as SKA, where many stable pulsars 
are expected to be discovered all over the sky. 

We begin by fixing the total coherent SNR at which the 
GW signal is observed , and we set SNR= 10, regard- 
less of the number of pulsars in the array, and explore 
the dependence of the results on the number of pulsars 
M in the range 3-to-lOOO. We then consider a fiducial 
'SKA-configuration' by fixing the total number of pul- 
sars to M = 100, and we explore how the results depend 
on the SNR for values 5 < SNR < 100. Throughout 
this analysis we assume that the timing noise is exactly 
the same for each pulsar and that the observations of 
each neutron star cover the same time span. The rela- 
tive contribution of each of the pulsars in the PTA to the 
SNR is therefore solely dictated by the geometry of the 
system pulsar-Earth-source, that is the specific value of 
the beam patter function (0, cj), ip). In total we con- 
sider 14 M-SNR combinations, and for each of them we 
generate 2.5 x lO^-to-1.6 x 10^ (depending on the total 
number of pulsars in the array) random sources in the 
sky. Each source is determined by the seven parameters 
described by Eq. pS]) . which, in all the Monte Carlo 
simulations presented from now on, are chosen as follow. 
The angles 9 and (f) are randomly sampled from a uni- 
form distribution in the sky; $o f^nd tp are drawn from a 
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FIG. 7: Median expected statistical error on the source pa- 
rameters. Each point (asterisk or square) is obtained by aver- 
aging over a large Monte Carlo sample of MBHBs (it ranges 
from 2.5 x 10'* when considering 1000 pulsars to 1.6 x 10^ when 
using 3 pulsars). In each panel, solid lines (squares) represent 
the median statistical error as a function of the total coherent 
SNR, assuming 100 randomly distributed pulsars in the sky; 
the thick dashed lines (asterisks) represent the median statis- 
tical error as a function of the number of pulsars M for a fixed 
total SNR= 10. In this latter case, thin dashed lines label the 
25"^ and the 75"^ percentile of the error distributions. 



uniform distribution over their relevant intervals, [0,27r] 
and [0,7r] respectively; t is sampled according to a proba- 
bility distribution = sint/2 in the interval [0,7r] and 
the frequency is fixed at / = 5 x 10~* Hz. Finally the 
amplitude R is set in such a way to normalise the signal 
to the pre-selected value of the SNR. For each source we 
generate M pulsars randomly located in the sky and we 
calculate the Fisher information matrix and its inverse as 
detailed in Section Hill We also performed trial runs con- 
sidering / = 10^^ Hz and / = 10^® Hz (not shown here) 
to further cross-check that the results do not depend on 
the actual GW frequency. 

Fig. [7] shows the median statistical errors as a function 
of M and SNR for all the six relevant source's parameters 
{6 and (p arc combined into the single quantity Af2, ac- 
cording to Eq. ([5^). Let us focus on the M dependence 
at a fixed SNR= 10. The crucial astrophysical quantity 
is the sky location accuracy, which ranges from « 3000 
deg^ for M = 3 ~ approximately 10% of the whole sky 
- to « 20 deg2 for M = 1000. A PTA of 100 pulsars 
would be able to locate a MBHB within a typical error 
box of K, 40 dcg^. The statistical errors for the other pa- 
rameters are very weakly dependent on M for M ^ 20. 
The fractional error in the source amplitude is typically 



FIG. 8: Distributions normalised to unity of the size of the 
error-box in the sky assuming an isotropic random distribu- 
tion of pulsars in the array. Upper panel: from right to left 
the number of pulsars considered is M = 3, 5, 20, 100, 1000, 
and we fixed a total SNR= 10 in all cases. Lower panel: from 
right to left we consider SNR= 5, 10, 20, 50, 100, and we fixed 
M = 100. 



« 30%, which unfortunately prevents to constrain an as- 
trophysically meaningful "slice" in the Ai — plane. 
The frequency of the source, which in this case was cho- 
sen to be / = 5 X 10^^ Hz, is determined at a ^ 0.1 
nHz level. Errors in the inclination and polarization an- 
gles arc typically « 0.3 rad, which may provide useful 
information about the orientation of the binary orbital 
plane. 

All the results have the expected scaling with respect 
to the SNR, i.e. Afl oc 1/SNR^, and for all the other 
parameters shown in Fig. [7] the uncertainties scale as 
1/SNR. A typical source with a SNR= 100 (which our 
current astrophysical understanding suggests to be fairly 
unlikely, see SVV) would be located in the sky within an 
error box ;$ 1 deg^ for M ^ 10, which would likely en- 
able the identification of any potential electro-magnetic 
counterpart. 

Distributions (normalised to unity) of AO are shown 
in Fig. m The lower panel shows dependence on SNR 
(at fixed number of pulsars in the PTA, here set to 100), 
whose effect is to shift the distributions to smaller val- 
ues of An as the SNR increases, without modifying the 
shape of the distribution. The upper panel shows the 
effectiveness of triangulation; by increasing the number 
of pulsars at fixed coherent SNR, not only the peak of 
the distribution shifts towards smaller values of AH, but 
the whole distribution becomes progressively narrower. If 
they yield the same SNR, PTAs containing a larger num- 
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FIG. 9: Median statistical error in the source's parameter 
estimation as a function of the sky-coverage of the pulsar dis- 
tribution composing the array. Each triangle is obtained av- 
eraging over a Monte Carlo generated sample of 1.6 x 10^ 
sources. In each panel, solid lines (triangles) represent the 
median error, assuming M = 100 and a total SNR= 10 in the 
array; thin dashed lines label the 25*'' and the TS'*" percentile 
in the statistical error distributions. 



ber of pulsars (sufficiently evenly distributed in the sky) 
with higher intrinsic noise are more powerful than PTAs 
containing fewer pulsars with very good timing stability, 
as they allow a more accurate parameter reconstruction 
(in particular for sky position) and they minimise the 
chance of GW sources to be located in "blind spots" in 
the sky (see next Section). 



C. Anisotropic distribution of pulsars 

The sky distribution of the pulsars in a PTA is not 
necessarily isotropic. This is in fact the case for present 
PTAs, and it is likely to remain the norm rather than the 
exception, until SKA comes on-line. It is therefore useful 
- as it also sheds new light on the ability of reconstruct- 
ing the source parameters based on the crucial location 
of the pulsars of the array with respect to a GW source ~ 
to explore the dependency of the results on what we call 
the "PTA sky coverage" A57pta, i-e. the minimum solid 
angle in the sky enclosing the whole population of the 
pulsars in the array. We consider as a study case a 'polar' 
distribution of 100 pulsars; the location in the sky of each 
pulsar is drawn from a uniform distribution in </> and cos 
with parameters in the range (j) G [0, 2tt] and 6 € [0, ^max] , 
respectively. We then generate a random population of 
GW sources in the sky and proceed exactly as we have de- 



scribed in the previous section. We consider six different 
values of AfipxA, progressively increasing the sky cover- 
age. We choose 0„iax = 7r/12, tt/6, 7r/4, tt/3, n/2, tt corre- 
sponding to ArjpTA = 0.21, 0.84, 1.84, tt, 27r, 47r srad. As 
we are interested in investigating the geometry effects, 
we fix in each case the total optimal SNR to 10. We ded- 
icate the next section to consider specifically the case of 
the 20 pulsars that arc currently part of the Parkcs PTA. 

The median statistic errors on the source parameters 
as a function of the PTA sky coverage are shown in Fig. 
|9l As one would expect, the errors decrease as the sky 
coverage increases, even if the SNR is kept constant. This 
is due to the fact that as the pulsars in the array populate 
more evenly the sky, they place increasingly more strin- 
gent constraints on the relative phase differences amongst 
the same GW signal measured at each pulsar, which de- 
pends on the geometrical factors F^'^ . The most im- 
portant effect is that the sky position is pinned down 
with greater accuracy; at the same time, correlations be- 
tween the sky location parameters and other parameters, 
in particular amplitude and inclination angle are reduced. 
An scales linearly (at fixed SNR) with A51pta, but the 
others parameters do not experience such a drastic im- 
provement. The statistical uncertainty on the amplitude 
improves as ^/AQp^A for AJIpta ^ 1 srad, then satu- 
rates. All the other parameters are much less sensitive 
to the sky coverage, showing only a mild improvement (a 
factor < 2) with increasing AfipxA up to '--^ 1 srad. 

When one considers an anisotropic distribution of pul- 
sars, the median values computed over a random uniform 
distribution of GW sources in the sky do not carry how- 
ever the full set of information. In particular the error- 
box in the sky strongly depends on the actual source 
location. To show and quantify this effect, we use the 
outputs of the Monte Carlo runs to build sky maps of 
the median of Aft that we shown in Fig. [TOl When 
the pulsars are clustered in a small AJIpta, the proper- 
ties of the signals coming from that spot in the sky (and 
from the diametrically opposite one) are more suscep- 
tible to small variations with the propagation direction 
(due to the structure of the response functions and 
F^); the sky location can then be determined with a 
much better accuracy, AJ7 ~ 2 deg^. Conversely, trian- 
gulation is much less effective for sources located at right 
angles with respect to the bulk of the pulsars. For a po- 
lar Ai7pTA = 0.21 srad, we find a typical Ail > 5000 
srad for equatorial sources; i.e., their sky location is ba- 
sically undetermined. Increasing the sky coverage of the 
array, obviously mitigates this effect, and in the limit 
AfipTA ~ 47r srad (which correspond to an isotropic pul- 
sar distribution), we find a smooth homogeneous skymap 
without any recognizable feature (bottom right panel of 
Fig. [TU)) . In this case the sky location accuracy is in- 
dependent on the source sky position and, for M = 100 
and SNR = 10 we find An ~ 40 deg^. Fig. [TT] shows 
the normalized distributions of the statistical errors cor- 
responding to the six skymaps shown in Fig. 1101 It is 
interesting to notice the bimodality of the distribution for 



15 









'^100 


— 100 




300 



200 



100 




■ 50 
45 
40 
35 



FIG. 10: Sky maps of the median sky location accuracy for an anisotropic distribution of pulsars in the array. Contour plots 
are generated by dividing the sky into 1600 (40 x 40) cells and considering all the randomly sampled sources falling within each 
cell; SNR= 10 is considered. The pulsar distribution progressively fills the sky starting from the top left, eventually reaching 
an isotropic distribution in the bottom right panel (in this case, no distinctive features are present in the sky map). In each 
panel, 100 black dots label an indicative distribution of 100 pulsars used to generate the maps, to highlight the sky coverage. 
Labels on the contours refer to the median sky location accuracy expressed in square degrees, and the color-scale is given by 
the bars located on the right of each map. 



intermediate values of AfJpTA; due to the fact that there 
is a sharp transition between sensitive and non sensitive 
areas in the sky (this is particularly evident looking at 
the contours in the bottom left panels of Fig. [TU]) . 



We also checked another anisotropic situation of po- 
tential interest: a distribution of pulsars clustered in the 
Galactic plane. We considered a distribution of pulsars 
covering a ring in the sky, with is randomly sampled in 
the interval [0,27r] and latitude in the range [— 7r/12, 7r/12] 
around the equatorial plane, corresponding to a solid an- 
gle of AfipTA = 3.26 srad. Assuming a source SNR= 10, 
the median statistical error in the source sky location 
is ~ 100 deg^, ranging from ~ 10 deg^ in the equato- 
rial plane, to ^ 400 deg^ at the poles. Median errors 
on the other parameters are basically the same as in the 
isotropic case. 



D. The Parkes Pulsar Timing Array 

We finally consider the case that is most relevant to 
present observations: the potential capabilities of the 
Parkes Pulsar Timing Array. The goal of the survey is to 
monitor 20 milli-second pulsars for five years with timing 
residuals « 100 ns [l[. This may be sufficient to enable 
the detection of the stochastic background generated by 
the whole population of MBHBs [28|, but according to 
our current astrophysical understanding (see SVV) it is 
unlikely to lead to the detection of radiation from indi- 
vidual resolvable MBHBs, although there is still a non- 
negligible chance of detection. It is therefore interesting 
to investigate the potential of such a survey. 

In our analysis we fix the location of the pulsars in 
the PTA to the coordinates of the 20 milli-second pul- 
sars in the Parkes PTA, obtained from [s^]; however for 
this exploratory analysis we set the noise spectral den- 
sity of the timing residuals to be the same for each pulsar. 
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FIG. 11: Normalized distributions of the statistical errors 
in sky position accuracy corresponding to the six sky maps 
shown in Fig. [10] Each distribution is generated using a 
random subsample of 2.5 x lO'' sources. 



i.e. we do not take into aecount the different timing sta- 
bility of the pulsars. We then generate a Monte Carlo 
sample of GW sources in the sky with the usual proce- 
dure. We consider two different approaches. Firstly, we 
explore the parameter estimation accuracy as a function 
of the GW source sky location for selected fixed array 
coherent SNRs (5, 10, 20, 50 and 100). Secondly, we fix 
the source chirp mass, frequency and distance (so that 
the sky and polarization averaged coherent SNR is 10) 
and we explore the parameter estimation accuracy as a 
function of the sky location. Skymaps of statistical er- 
ror in the sky location are shown in Fig. [T^l In the top 
panel we fix the SNR= 10, independently on the source 
position in the sky; the median error in the sky location 
accuracy is Afi ^ 130 deg^, but it ranges from ^ 10 deg^ 
to ~ 400 deg^ depending on the source's sky location. 
The median statistical errors that affect the determina- 
tion of all the other source parameters are very similar 
to those for the isotropic pulsar distribution case when 
considering M = 20, since the pulsar array covers al- 
most half of the sky, see Fig. [HI In the bottom panel, 
we show the results when we fix the source parameters, 
and therefore, the total SNR in the array does depend 
on the source sky location. In the southern hemisphere, 
where almost all the pulsars are concentrated, the SNR 
can be as high as 15, while in the northern hemisphere it 
can easily go below 6. The general shape of the skymap 
is mildly affected, and shows an even larger imbalance 
between the two hemispheres. In this case, the median 
error is Af2 ^ 160 deg^, ranging from ~ 3 deg^ to ~ 900 
deg^. It is fairly clear that adding a small ( ^ 10) number 



pulsars in the northern hemisphere to the pulsars already 
part of the Parkes PTA would significantly improve the 
uniformity of the array sensitivity and parameter estima- 
tion capability, reducing the risk of potentially detectable 
GW sources ending up in a "blind spot" of the array. 



V. CONCLUSIONS 

In this paper we have studied the expected uncertain- 
ties in the measurements of the parameters of a mas- 
sive black hole binary systems by means of gravitational 
wave observations with Pulsar Timing Arrays. We have 
investigated how the results vary as a function of the 
signal-to-noise ratio, the number of pulsars in the array 
and their location in the sky with respect to a gravita- 
tional wave source. Our analysis is focused on MBHBs in 
circular orbit with negligible frequency evolution during 
the observation time ("monochromatic sources"), which 
we have shown to represent the majority of the observ- 
able sample, for sensible models of sub-parsec MBHB 
eccentricity evolution. The statistical errors are evalu- 
ated by computing the variance- covariance matrix of the 
observable parameters, assuming a coherent analysis of 
the Earth-terms only produced by the timing residuals 
of the pulsars in the array (see Section II B). 

For a fiducial case of an array of 100 pulsars randomly 
distributed in the sky, assuming a coherent total SNR 
= 10, we find a typical error box in the sky Ail « 40 
deg^ and a fractional amplitude error of « 0.3. The lat- 
ter places only very weak constraints on the chirp mass- 
distance combination M^^^ /Dl. At fixed SNR, the typ- 
ical parameter accuracy is a very steep function of the 
number of pulsars in the PTA up to ~ 20. For PTAs 
containing more pulsars, the actual gain becomes pro- 
gressively smaller because the pulsars "fill the sky" and 
the effectiveness of further triangulation weakens. We 
also explored the impact of having an anisotropic dis- 
tribution of pulsars finding that the typical source sky 
location accuracy improves linearly with the array sky 
coverage. For the specific case of the Parkes PTA where 
all the pulsars are located in the southern sky, the sensi- 
tivity and sky localisation are significantly better (by an 
order of magnitude) in the southern hemisphere, where 
the error-box is < 10 deg^ for a total coherent SNR = 
10. In the northern hemisphere, the lack of monitored 
pulsars prevent a source location to be in an uncertainty 
region < 200 deg^ . The monitoring of a handful of pul- 
sars in the northern hemisphere would significantly in- 
crease both the SNR and the parameter recovery of GW 
sources, and the International PTA [3| will provide such 
a capability in the short term future. 

The main focus of our analysis is on the sky locali- 
sation because sufficiently small error-boxes in the sky 
may allow the identification of an electro-magnetic coun- 
terpart to a GW source. Even for error-boxes of the order 
of tens-to-hundreds of square degrees (much larger than 
e.g. the typical LISA error-boxes |53l - [55| ). the typical 
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FIG. 12: Sky maps of the median sky location accuracy for the Parkes PTA. Contour plots are generated as in Fig. 1101 Top 
panel: we fix the source SNR= 10 over the whole sky; in this case the sky position accuracy depends only on the different 
triangulation effectiveness as a function of the source sky location. Bottom panel: we fix the source chirp mass and distance to 
give a sky and polarization averaged SNR= 10, and we consistently compute the mean SNR as a function of the sky position. 
The sky map is the result of the combination of triangulation efficiency and SNR as a function of the source sky location. The 
color-scale is given by the bars on the right, with solid angles expressed in deg^. 



sources are expected to be massive {M ^ IO^Mq) and 
at low redshift (z ^ 1.5), and therefore the number of 
associated massive galaxies in the error-box should be 
Hmited to a few hundreds. Signs of a recent merger, like 
the presence of tidal tails or irregularities in the galaxy 
luminosity profile, may help in the identification of po- 
tential counterparts. Furthermore, if nuclear activity is 
present, e.g. in form of some accretion mechanism, the 
number of candidate counterparts would shrink to an 
handful, and periodic variability [56j could help in as- 
sociating the correct galaxy host. We arc currently in- 
vestigating the astrophysical scenarios and possible ob- 
servational signatures, and we plan to come back to this 
important point in the future. The advantage of a coun- 
terpart is obvious: the redshift measurement would allow 
us, by assuming the standard concordance cosmology, to 
measure the luminosity distance to the GW source, which 
in turn would break the degeneracy in the amplitude of 
the timing residuals R cx M^^^ /{DlJ^^^) between the 
chirp mass and the distance, providing therefore a direct 



measure of ^A. 

The study presented in this paper deals with 
monochromatic signals. However, the detection of MB- 
HBs which exhibit a measurable frequency drift would 
give significant payoffs, as it would allow to break the 
degeneracy between distance and chirp mass, and en- 
able the direct measurement of both parameters. Such 
systems may be observable with the Square-Kilometre- 
Array. In the future, it is therefore important to extend 
the present analysis to these more general signals. How- 
ever, as the frequency derivative has only modest cor- 
relations with the sky position parameters, we expect 
that the results for the determination of the error-box 
in the sky discussed in this paper will still hold. A fur- 
ther extension to the work is to consider MBHBs charac- 
terised by non-negligible eccentricity, which is currently 
in progress. Another extension to our present study is to 
consider both the Earth- and pulsar-terms in the analysis 
of the data and the investigation of the possible benefits 
of such scheme, assuming that the pulsar distance is not 
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known to sufficient accuracy. This also raises the issue 
of possible observation campaigns that could yield an ac- 
curate (to better than 1 pc) determination of the pulsar 
distances used in PTAs. In this case the use of the pulsar- 
term in the analysis would not require the introduction 
of (many more) unknown parameters and would have the 
great benefit of breaking the degeneracy between chirp 
mass and distance. 

The final world of caution goes to the interpretation 
of the results that we have presented in the paper. The 
approach based on the computation of the Fisher Infor- 
mation matrix is powerful and straightforward, and is 
justified at this stage to understand the broad capabili- 
ties of PTAs and to explore the impact on astronomy of 
different observational strategies. However, the statisti- 
cal errors that we compute are strictly lower limits to the 
actual errors obtained in a real analysis; the fact that at 
least until SKA comes on line, a detection of a MBHB 
will be at a moderate-to-low SNR should induce caution 
in the way in which the results presented here are inter- 
preted. Moreover, in our current investigation, we have 
not dealt with a number of important effects that in real 
life play a significant role, such as different calibrations 
of different data sets, the change of systematic factors 
that affect the noise, possible non-Gaussianity and non- 
stationarity of the noise, etc. These (and other) impor- 
tant issues for the study of MBHBs with PTAs should 
be addressed more thoroughly in the future by perform- 
ing actual mock analyses and developing suitable analysis 
algorithms. 



shift at time t, Eq. ([3|) is therefore 



Appendix A: The pulsar timing array response to 
gravitational waves 

In this Appendix we derive the PTA response to a GW 
signal from a deterministic source characterised by a met- 
ric perturbation 



hit,xo) = A+{t,xo)F+{Cl) - iA^{t,xo)F'' (il) 



(Al) 

The observed signal at the PTA output generated by 
the fractional frequency shift comes from the contribu- 
tion at the "Earth" (or more precisely the solar system 
barycentre, SSB) and at the pulsar. We therefore need 
to compute the metric perturbation at (tssb,a;ssb) and 
{t-p^Xp). For simplicity (but without loss of generality), 
let us make a choice of coordinates such that: 



^ *^ssb — : 



t 



where 



(A3) 



is the time delay between the pulsar and the Earth terms, 
given by the light-travel-time. The relative frequency 



A+{t - t)F+ (Cl) - iA^{t - t)F'' (Cl) 



A+{t)F+{n)-iAy,{t)F''{n) 



(A4) 



As one expects, the response to a GW is identically zero 
for GWs propagating along the Earth-pulsar direction, 
when Cl = zLp. In fact, when GWs propagate in the 
direction opposite to the radio signal, Q, ■ p = 1, we have 
rh-p = h-p = and as a consequence F^{Ct) — F^ {Ct) = 
0. This is due to the transverse nature of GWs and there 
is no frequency shift of the electro-magnetic wave. On 
the other hand, when GWs propagate along the same 
direction of the electromagnetic waves, Cl ■ p = —1 and 
F+(f2) = cos(2'i/') and F^(f2) = sm{2ip) are non-zero, 
finite functions of the polarization angle ■0. In this case 
z{t, Cl) is still zero because r = and the Earth and 
pulsar terms are identical and cancel out. This is known 
as the surfing effect, as the GWs surf with the electro- 
magnetic waves. 

We have shown in Section |lT] that GWs emitted 
by MBHBs and observable with PTA will be quasi- 
monochromatic signals slowly drifting in time. More 
specifically, for the astrophysical population that one ex- 
pects to detect (see SVV), the frequency derivative sat- 
isfies (for the vast majority of the signals) the condition 



/o 



(A5) 



The consequence of the inequality above is that the tim- 
ing residuals observed at the PTA output from any pul- 
sar a will consist of two quasi-monochromatic signals at 
different, but essentially constant frequencies during the 
observation time f(t) and f(t — Tq,), where we have ex- 
plicitly emphasised with Tq the fact that the pulsar-term 
differs from pulsar to pulsar, and will be sitting at a dif- 
ferent frequency determined by the source-Earth-pulsar 
relative angle Q-pa and the poorly constrained pulsar dis- 
tance La- Given then the fact that the noise dominates 
the signal, <C n^, only adding the contributions to all 
the pulsars will provide enough SNR, and one can there- 
fore concentrate the analysis only on the Earth-terms, all 
of which depend only on the 7 parameters of the source, 
and ignore the contribution from the pulsar terms. The 
latter in fact depend on the M unknown parameters L^s, 
which determine Tq, and therefore phase, frequency and 
amplitude of the pulsar term. If L^s were known - an 
uncertainty smaller than 1 rad on the phase contribution 
requires the distance to a pulsar to be known to better 
(A2) than ss 0.1 (//lOnHz)"^ (1 + Cl-p)~^ pc - one can coher- 
ently lock all the phases of the pulsar terms increasing 
the total SNR of the detection. This would also pro- 
vide a direct measurement of /, allowing to break the 
distance-chirp mass degeneracy. 

The pulsar-term however may conjure to cancel the 
Earth-term for specific source-Earth-pulsar angles. To 
have an order of magnitude estimate of this effect, let 
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us conservatively assume that the Earth and pulsar term 
can be fully resolved if their frequency separation is larger 
than one frequency resolution bin of width 1 /T (we have 
actually shown in Section IIVI that the frequency can be 
resolved with sub- frequency resolution accuracy). As- 
suming (again conservatively) a linear frequency shift due 
to radiation reaction, this condition can be expressed as 



i 1 
Jot > ^ 



(A6) 



Substituting Eq. (|A3|) in the previous expression, the pre- 
vious condition can be expressed as 



(A7) 



The cancellation takes place when f2 • p ~ 1 and we can 
therefore approximate 1 + ■ p as 



1 + ■ p = 1 + cos(7r 



1 



(A8) 

Using a typical pulsar distance of 1 kpc and a GW fre- 



quency of 10 nHz we obtain: 



10nHz\ 



1/2 



1 kpc 
L 



1/2 



deg. 



(A9) 



This means that whenever the separation on the sky of 
a GW source and a pulsar is larger than a few degrees, 
then ignoring the contribution from the pulsar term does 
not affect the analysis. With existing PTAs, this is 
surely a safe assumption unless one encounters a very 
unlucky case. Nonetheless, we keep this into account in 
the Monte Carlo simulations whose results are reported 
in Section ITVl Wc choose 50nHz for the source frequency. 
Then, for every PTA array configuration, we place all 
the pulsars either at 1 kpc or at 5 kpc (in Section ITVl we 
show results for La = 5 kpc, but the results for La = 1 
kpc are basically identical); while for the Parkes PTA 
we took each individual pulsar distance from the ATNF 
catalogue [52[. Finally, for each pulsar we compute the 
pulsar-Earth-source angle and if the resulting value is 
smaller than the one obtained from Eq. (jA9p we discard 
that particular pulsar contribution from the analysis; as 
66 is very small, the impact on the analysis is minimal. 
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